rm(list=ls())  # Clear the workspace

# Create a log file
sink(file="Avery_Brevoort_Final_Graph.log",
     append=FALSE,
     type="output",
     split=TRUE)

# Read in the data
load("REStatSubmission_FRB.RData")

# Only retain the observations that will be displayed (Census Tracts with relative incomes
# of between 60 and 140 percent)
myData <- subset(myData, (relinc>60 & relinc<=140) )

# Create a smaller data frame that contains only the variables that will be plotted
graphData <- data.frame( count=myData$COUNT,
                         delrate=myData$delrate,
                         msa=myData$MSA,
                         relinc=floor(myData$relinc),
                         hidti_20046=myData$hidti_20046,
                         hirate_20046=myData$hirate_20046,
                         gse_20046=myData$solf_20046,
                         shr.dep.in_20046=myData$shr.dep.in_20046,
                         pur.dep.in_20046=myData$pur.dep.in_20046,
                         loans_20046=myData$loans_20046
                         )

#######################################################################
###  Figure 1 (This is the only figure in the paper)
#######################################################################

# Set up the output file
png("DiscontinuityGraphs.png",
    width=6.5,
    height=8,
    units="in",
    res=600)

# Set up display as a matrix of 6 graphs (3 x 2) and set margins around the graphs
par(mfrow=c(3,2),
    mar=c(4,3,2,1) )

# This function calculates weighted means by specified group.  The groups used in the 
# graphs will be the levels of relative income between 60 and 140.
ave.mean <- function( x, w, group ) {
  numerator <- aggregate( x*w, by=group, FUN=sum )
  denominator <- aggregate( w, by=group, FUN=sum )
  group.id <- numerator[,1]
  group.mean <- numerator[,2]/denominator[,2]
  return( cbind(group.id, group.mean ) )
}

# This function demeans a variable (x), using a weighted mean (with weight equal to w),
# by group.  The group used in the graphs will be the relative income levels between 60 
# and 140 percent.  Adds in the grand mean for the variable.
demean.variable <- function( x, w, group ) {
  grand.mean <- weighted.mean( x=x, w=w )
  group.mean <- ave.mean( x, w=w, group=list(group) )
  dimnames(group.mean)[[2]] <- c("group","group.mean")
  group.mean <- data.frame( group.mean )
  tempData <- data.frame( x=x, w=w, group=group )
  tempData <- merge( tempData, group.mean, by="group" )
  tempData$grand.mean <- grand.mean
  retVar <- tempData$x - tempData$group.mean + tempData$grand.mean
  return(retVar)
}

# Demean the six variables that will be plotted
graphData$delrate.dm <- demean.variable( x=graphData$delrate,
                                         w=graphData$count,
                                         group=graphData$msa )
graphData$hidti_20046.dm <- demean.variable( x=graphData$hidti_20046,
                                             w=graphData$loans_20046,
                                             group=graphData$msa )
graphData$hirate_20046.dm <- demean.variable( x=graphData$hirate_20046,
                                              w=graphData$loans_20046,
                                              group=graphData$msa )
graphData$gse_20046.dm <- demean.variable( x=graphData$gse_20046,
                                           w=graphData$loans_20046,
                                           group=graphData$msa )
graphData$shr.dep.in_20046.dm <- demean.variable( x=graphData$shr.dep.in_20046,
                                                  w=graphData$loans_20046,
                                                  group=graphData$msa )
graphData$pur.dep.in_20046.dm <- demean.variable( x=graphData$pur.dep.in_20046,
                                                  w=graphData$loans_20046,
                                                  group=graphData$msa )

#############################################################
### PANEL A - Delinquency Rate
#############################################################

# Calculate the average of the demeaned variable being plotted 
# at each relative income level between 60 and 140.
temp <- ave.mean( x=graphData$delrate.dm,
                  w=graphData$count, 
                  group=list(graphData$relinc) )
temp <- subset( data.frame(temp), 
                group.id>60 & group.id<=140 )

plot(x=temp[,1],
     y=temp[,2],
     xlim=c(60,140),
     #ylim=c( min(temp[,2])*0.9, max(temp[,2]*1.1) ),
     type="o",
     pch=20,
     main="(a) Delinquency Rate",
     bty="l",
     ylab="",
     xaxt="n",
     yaxt="n",
     xlab="")
legend("topright",
       c("CRA Threshold","GSE Threshold"),
       lty=c(3,2),
       cex=1,
       xjust=1,
       yjust=0.5)
axis(2,las=1)
axis(1,las=1,at=c(60,70,80,90,100,110,120,130,140))
abline(v=90,lty=2)
abline(v=80,lty=3)

#############################################################
### PANEL B - High Payment-to-Income Share
#############################################################

# Calculate the average of the demeaned variable being plotted 
# at each relative income level between 60 and 140.
temp <- ave.mean( x=graphData$hidti_20046.dm,
                  w=graphData$loans_20046, 
                  group=list(graphData$relinc) )
temp <- subset( data.frame(temp), 
                group.id>60 & group.id<=140 )

plot(x=temp[,1],
     y=temp[,2],
     xlim=c(60,140),
     type="o",
     pch=20,
     main="(b) % High PTI",
     bty="l",
     ylab="",
     xaxt="n",
     yaxt="n",
     xlab="")
axis(2,las=1)
axis(1,las=1,at=c(60,70,80,90,100,110,120,130,140))
abline(v=90,lty=2)
abline(v=80,lty=3)

#############################################################
### PANEL C - Share of Higher-Priced Loans
#############################################################

# Calculate the average of the demeaned variable being plotted 
# at each relative income level between 60 and 140.
temp <- ave.mean( x=graphData$hirate_20046.dm,
                  w=graphData$loans_20046, 
                  group=list(graphData$relinc) )
temp <- subset( data.frame(temp), 
                group.id>60 & group.id<=140 )

plot(x=temp[,1],
     y=temp[,2],
     xlim=c(60,140),
     type="o",
     pch=20,
     main="(c) % Higher Priced",
     bty="l",
     ylab="",
     xaxt="n",
     yaxt="n",
     xlab="")
axis(2,las=1)
axis(1,las=1,at=c(60,70,80,90,100,110,120,130,140))
abline(v=90,lty=2)
abline(v=80,lty=3)

#############################################################
### PANEL D - GSE Sales (2004-2006)
#############################################################

# Calculate the average of the demeaned variable being plotted 
# at each relative income level between 60 and 140.
temp <- ave.mean( x=graphData$gse_20046.dm,
                  w=graphData$loans_20046, 
                  group=list(graphData$relinc) )
temp <- subset( data.frame(temp), 
                group.id>60 & group.id<=140 )

plot(x=temp[,1],
     y=temp[,2],
     xlim=c(60,140),
     type="o",
     pch=20,
     main="(d) GSE Sales",
     bty="l",
     ylab="",
     xaxt="n",
     yaxt="n",
     xlab="")
axis(2,las=1)
axis(1,las=1,at=c(60,70,80,90,100,110,120,130,140))
abline(v=90,lty=2)
abline(v=80,lty=3)

###########################################################################
### PANEL E - Lending by Depositiories Within Assessment Areas (2004-2006)
###########################################################################

# Calculate the average of the demeaned variable being plotted 
# at each relative income level between 60 and 140.
temp <- ave.mean( x=graphData$shr.dep.in_20046.dm,
                  w=graphData$loans_20046, 
                  group=list(graphData$relinc) )
temp <- subset( data.frame(temp), 
                group.id>60 & group.id<=140 )

plot(x=temp[,1],
     y=temp[,2],
     xlim=c(60,140),
     type="o",
     pch=20,
     main="(e) Lending: Depository - In",
     bty="l",
     ylab="",
     xaxt="n",
     yaxt="n",
     xlab="Relative Tract Income")
axis(2,las=1)
axis(1,las=1,at=c(60,70,80,90,100,110,120,130,140))
abline(v=90,lty=2)
abline(v=80,lty=3)

#####################################################################################
### PANEL F - Share of Purchases by Depositories within Assessment Areas (2004-2006)
#####################################################################################

# Calculate the average of the demeaned variable being plotted 
# at each relative income level between 60 and 140.
temp <- ave.mean( x=graphData$pur.dep.in_20046.dm,
                  w=graphData$loans_20046, 
                  group=list(graphData$relinc) )
temp <- subset( data.frame(temp), 
                group.id>60 & group.id<=140 )

plot(x=temp[,1],
     y=temp[,2],
     xlim=c(60,140),
     type="o",
     pch=20,
     main="(f) Purchases: Depository - In",
     bty="l",
     ylab="",
     xaxt="n",
     yaxt="n",
     xlab="Relative Tract Income")
axis(2,las=1)
axis(1,las=1,at=c(60,70,80,90,100,110,120,130,140))
abline(v=90,lty=2)
abline(v=80,lty=3)

dev.off()

sink()  # Close log file

